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1 INTRODUCTION 



The field of solar system dynamics, fias a timeline of dis- 
coveries th at is related to the computational power avail- 
able (e.g., iMorbideiiil l200ll ). As computational power has 
increased over time, so to has our ability to more accu- 
rately simulate more complex systems. The interaction be- 
tween planets and planetesimals subsequent to formation is 
a well posed n-body problem but displays remarkable com- 
plexity even in the absence of collisions, including resonance 
capture, planetary migration, and heating and scattering of 
planetesimals. In this paper we attempt to more accurately 
simulate this system by using the increased computational 
power and on board memory recently available on video 
graphics cards. 

Planet-planetesi mal scattering can lea d to planet 

migra tio n (e.g., Fernandez fc id 1984 1 Malhotra 
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ABSTRACT 

We simulate planet migration caused by interactions between planets and a planetes- 
imal disk. We use an N-body integrator optimized for near-Keplerian motion that 
runs in parallel on a video graphics card, and that computes all pair-wise gravita- 
tional interactions. We find that the fraction of planetesimals found in mean motion 
resonances is reduced and planetary migration rates are on average about 50% slower 
when gravitational interactions between the planetesimals are computed than when 
planetesimal self-gravity is neglected. This is likely due to gravitational stirring of the 
planetesimal disk that is not present when self-gravity is neglected that reduces their 
capture efficiency because of the increased particle eccentricity dispersion. We find 
that migration is more stochastic when the disk is self-gravitating or comprised of 
more massive bodies. Previous studies have found that if the planetesimal disk den- 
sity is below a critical level, migration is "damped" and the migration rate decays 
exponentially, otherwise it is "forced" and the planet's migration rate could accelerate 
exponentially. Migration rates measured from our undamped simulations suggest that 
the migration rate saturates at a level proportional to disk density and subsequently 
is approximately power law in form with time. 

planetesimals. However, because of the 0{N^) required 
computational intensity previous simulations have neces- 
sarily neglected the gravitational interaction s between the 
planetesimals (e.g., 'Hahn fc Malhotral 19991: iGomes et al.l 
2004; Hahn & Malhotra 2005; Thom mes et alj 120081 ). This 
means that collective gravitational effects and gravitational 
self stirring in the disk have been ignored. Here we take all 
interparticle forces into account to explore what differences 
might be seen in simulations of migrating planets with and 
without planetesimal interactions. 



Levison fc Morbidelli 



19951: iHahn fc Malhotral Il999l: llda et al 



20031: iG omes et al.l |2004 



20051 : 



Hahn fc M alhotra 2005|; I Levison et al. 20081 ). Recent work 
on dusty circumstellar disks with clearings, suggest that 
plane ts are required to account for the morphology of the 
disk (|QuiUenll2006al : IWvattll2006l ). Th e planetesimal masses 
suggeste d by coUisional models (e.g.. 



[Pominik fc Decinll20oA iQuiUen et all 



Wvatt fc Dentll2002l: 



2007h and proximity 



of the hypothetical planet to the disks for systems such 
as Fomalhaut intimate that planetary migration could be 
taking place. To simulate planet migration due to planetes- 
imal scattering, planets must gravitationally interact with 



2 NUMERICAL SIMULATIONS 

The migration simulations were carried out using an N-body 
code running on a 4 node cluster running the "Mars Hill" 
Rocks 4.3 Linux operating system (a CentOS based distri- 
bution) that hosts 4 NVIDIA GeForce 8800 GTX graphics 
cards. The code was run on the graphics processing units 
(GPUs) residing on the graphics cards and is written with 
NVIDIA's CUDA (Compute Unified Device Architecture), 
a C-language development environment for CUDA enabled 
GPUs. CUDA is a GPGPU (General-purpose computing on 
graphics processing units) technology that allows a program- 
mer to use the C programming language to code algorithms 
for execution on the graphics processing unit. Its intention is 
to expose the hardware to the developer through a memory 
management model that encourages both constant stream- 
ing of data as well as massive parallelization. 
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2.1 A second order democratic heliocentric 

method symplectic integrator for the GPU 

We have mod ified the second ord er symplectic integrator in- 
troduced by jPuncan et al.|[T998l ). also known as the demo- 
cratic heliocentric method, so that it runs on a GPU. We 
have chosen the democratic heliocentric method because the 
force from the central body is separated from the integra- 
tion of the remainder of the particles and the coordinates do 
not depend on the order of the particles. This increases the 
accuracy of the integrator when there are large mass differ- 
ences and is particularly desirable when forced to work in 
floating point single precision. When we obtain video cards 
that can compute in double precision the code will improve 
in precision. The move to double precision will also allow us 
to extend the range of the planetesimal masses simulated. 

In heliocentric coordinates an d barycentric momenta 
l| Wisdom. Holman fc Toumalliggd ) the Hamiltonian of the 
system can be written 



H — Hsmi + Hkbp + Hint 



where 



2mo 



(1) 



(2) 



is a linear drift term and Pi are the barycentric momenta. 
Here mo is the central particle mass. The second term Hkbp 
is the sum of Keplerian Hamiltonians for all particles with 
respect to the central body, 

i — 

where Qi are the heliocentric coordinates and are conjugate 
to the barycentric momenta. Here rrn is the mass of the i-th 
particle and G is the gravitational constant. The interac- 
tion term contains all gravitational interaction terms except 
those to the central body. 
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The second order single timestep integrator advances 
with timestep t using evolution operators (e.g., lYoshidal 
1 19901 ) 

Esun Ekc'p Eint (t) Enep Esun (5) 

where we have reversed the order of the Keplerian evolu- 
tion and the interaction steps compared to that discussed 
by (jPuncan et al.|[l998l ). We have done this to reduce the 
total number of computations per timestep. The Keplerian 
advance requires 0{N) computations but interaction term 
requires 0{N^) computations. 

The drift evolution operator requires computation of 
the sum of the momenta. We have implemented this using a 
parallel reduction sum algorithm available with the NVIDIA 
CUD A Software Development Kit (SDK) 1.1 that is similar 
to the scan prefix sum algorithm ([Harris et al.ll2008h . 

The Keplerian step was implemented with / and g 
functions using the univer sal differential Kepler's equation 
l|Prussing fc Conwavl[l993l ) so that bound and unbound par- 
ticles can both be integrated with the same routine. The 



Keplerian evolution step is also done on the GPU with each 
thread computing the evolution for a separate particle. The 
dominant source of error is in the Keplerian evolution step 
and is due to the single floating point precision. These errors 
can cause a systematic radial drift that does not average to 
zero. To minimize errors caused by the single precision com- 
putation during the Kepler advances we chose / and g func- 
tions that maintain angular momentum conservation across 
each evolution step. The positions and velocities at a later 
time can be written in terms of those at an earlier time 



Xl 
Vl 



/xo + gvo 
fxo + gvo- 



(6) 



The angular momentum at the later time, Li can be written 
in terms of that at the earlier time, Lo, 



Li = {fa - 3/)(xo X vo) = {fg - gf)U- 



(7) 



Conservation of angular momentum yields the condition 

fi) - 9f = 1- (8) 

We utilize this formula to solve for one of the 4 functions 
reducing the inward radial drift resulting by the single pre- 
cision computation during the Keplerian advances. 

The interaction terms are computed on the GPU 
with all A*'^ force pairs evaluated explicitly in parallel. 
The algorithm is based on t he algorithm described by 
dNvland. Harris, fc Prin3l2008h . This algorithm takes ad- 
vantage of fast shared memory on board the GPU to si- 
multaneously compute all forces in a p x p tile of particle 
positions, where p is the number of threads chosen for the 
computation (typically 256). The total energy was evalu- 
ated with a kernel explicitly evaluating all A'^'^ pair potential 
energy terms, similar to that calculating all forces. 

After the change to heliocentric/barycentric coordi- 
nates, the position of the first coordinate corresponds to 
the center of mass and center of momentum. The trajectory 
of this particle need not be integrated. However it is con- 
venient to calculate the energy using all pair interactions 
including the central mass. The interaction term in the Ke- 
plerian part of the Hamiltonian can be computed at the 
same times as Hi„t if Qo is set to zero. Consequently we set 
Qo = Po = at the beginning of the computation. This is 
equivalent to working in the center of mass and momentum 
reference frame. Because we would like to be able to quickly 
check the total energy, we have chosen to keep the first par- 
ticle corresponding to the center of mass and momentum as 
the first element in the position and velocity arrays. During 
computation of Hint we set mo to zero so that force terms 
from the first particle are not computed. These are already 
taken into account in the evolution term corresponding to 
Hkep- The mass is restored during the energy sum com- 
putation as all potential energy terms must be calculated 
explicitly. 

The location of memory should be considered when run- 
ning routines as data transfer between the CPU and GPU 
can slow the code. The maximum theoretical throughput of 
PCI-express 16x bus technology, the interlink between the 
CPU and GPU on Intel processor based motherboards, is 
4 GB/s with a significantly higher latency than on device 
memory transfers. Depending on GPU design, theoretical 
memory transfer throughputs can approach 100 GB/s, as is 
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the case for the 8800 GTX which has a theoretical through- 
put of 86.4 GB/s. Reducing the number of data transfers 
between the CPU and GPU not only reduces the CPU clock 
time required for a calculation, but also reduces the number 
of calls which have higher latency. CPU and GPU global 
memory access is a significant bottleneck compared to the 
actual number of clock cycles required for a specific calcula- 
tion with the former containing a much higher performance 
penalty due to latency. For this reason, positions and veloc- 
ities for all particles are kept in global memory on board the 
GPU device. Arrays for both current and previous time step 
positions are required during the interaction step computa- 
tion. The particle positions and velocities are only trans- 
ferred back into host or CPU accessible memory to output 
data files. An additional float vector of length equal to the 
number of particles is allocated in global memory on the de- 
vice to compute the momentum sums used in the drift step 
computation. Shared memory on the GPU is used during 
the all pairs interaction computation step and during the 
reduction sum. Shared memory on the GPU drastically in- 
creases the speed of calculations because the use of shared 
memory has almost no latency penalty compared to that of 
both the CPU and the GPU global memory. By streaming 
information from global memory to shared memory, we are 
able to hide the latency of global memory, further increas- 
ing the computation speed. Though the maximum number 
of threads on the video card we used is 512, we found that 
register space on the GPU limited the interaction step com- 
putation to 256 threads per block due to its complexity. A 
more detailed review of NVIDIA GPU hardware and pro- 
gramming techniques can be found in the CUDA Program- 
ming Guide. 

We work in a lengthscale in units of the outermost 
planet's initial semi-major axis and with a timescale such 
that GM* — 1 where M, is the mass of the central star. 

2.2 Code Checks 

The following self-consistency checks on the code were per- 
formed. We checked that two body dynamics is preserved 
for all particles when all masses but the central one are 
zero. In this case the orbital elements (excepting the mean 
anomaly) are conserved at the level of the precision of com- 
putation. We checked that the energy computation com- 
puted in heliocentric/barycentric coordinates is consistent 
with that computed in Cartesian coordinates. We checked 
that the conversion to heliocentric/barycentric coordinates 
followed by inverse conversion returns the initial coordi- 
nates. We checked that the force and potential energy terms 
computed on the GPU are consistent with those computed 
on CPU. We checked that the N-body simulation without 
Keplerian integration conserves energy. We note that the 
order of some of the arguments in the SDKl.l nbody ker- 
nel were inconsistent with the correct order described by 
l|Nvland. Harris, fc Prinj|2008l ). To conserve energy we cor- 
rected the order of the arguments in the force computation 
step from that in the SDKl.l nbody kernel. We checked that 
integrating a two body system conserves energy. We ran a 
test case identical to that run bv lBromlev fc KenvonI (120061 ) 
and shown in their Figure 5 with similar results (though see 
discussion on energy conservation in the next section). This 
corresponds to 2 earth mass planets embedded in a uniform 



disk of 800 objects 1/100 of th e mass of the p l anets . This 
test case was also integrated bv lKokubo fc Ida (|l995l ). 



2.3 Simulation Checks 

The second set of checks was done by running full simula- 
tions in order to determine the smoothing length and time 
step requirements and restrictions. In one test we varied the 
time step significantly in order to measure the impact of 
larger time steps on our accuracy. A larger timestep would 
allow us to run the simulation for more orbital periods. We 
ran full lO** particle simulations over thousands of orbits 
and noted that varying the time step resulted in little dif- 
ference in the total energy measured. In our second test, 
we varied the softening length in 2 different sets of simula- 
tions. At extremely small softening length values we note 
that the simulated system becomes extremely stochastic. 
The reverse effect is exhibited when we make the smoothing 
length extremely larg e. This effect was previously noted in 
(|Kokubo fc Idalll995l ) and happens because our simulation 
does not adjust the time step to take into account collisions 
and close encounters. Since the simulation steps are small 
enough to allow "smooth" or nearly continuous updates, 
there is the possibility that a close encounter between two 
objects will be forced by the step which may have never oc- 
curred physically. If the softening length is made too large, 
the forces between objects that are having close encounters 
is reduced unrealistically. The sof tening length is chosen to 
help prevent this from happening. iKokubo fc Idal (|l995l ) se- 
lects a softening parameter set as 



where vh is the Hill radius of a planetesimal. This leads to 
a Hill radius which is roughly equivalent to the physical ra- 
dius of a protoplanet at lAU. Given these parameters, we 
chose a variety of smoothing lengths that were on the same 
order of magnitude and did not note a significant any signif- 
icant deviations between the results of the simulations. The 
smoothing length and timestep were not adjusted during the 
simulation. 



2.4 Benchmarks 

We measured the time to run our code on on the NVIDIA 
GeForce 8800 GTX graphics cards. The timing included the 
delay to send the job out to the individual nodes of the clus- 
ter, initial disk generation and input and output of files. We 
ran a simulation of a disk with 2 planets, one at semi-major 
axis a = 1 with the other interior to this. The planet masses 
were 10"'^ of that of the star. The planets were surrounded 
by a disk of particles with mass 10~^ of that of the central 
star. Time units were set such that GM* = 1. The time 
step was 0.1 or 0.1/(27r) = 0.0159 of the outer planet's or- 
bit. For 10* particles 1000 timesteps ran with 0.137 s per 
timestep. For 10^ particles 100 timesteps ran with 57s per 
timestep and similar error. The energy fraction error was 
10"'^ throughout the entire simulation and drifts downward 
systematically with time such that the system becomes in- 
creasingly bound. The systematic loss of energy and associ- 
ated radial inward drift we have determined is caused by the 
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single point floating precision for calculations done during 
the Kepler advances. 

2.5 Initial conditions 

In this paper we primarily discuss simulations with 2 planets 
and an exterior annulus of planetesimals. The initial plan- 
etesimal disk has a flat distribution in semi-major axis so 
that the surface density is proportional to 1/r for radius r. 
An isotropic disk is generated with initial eccentricity and 
inclination dispersion such that (e^) — 2{i^) and initial in- 
clination dispersion set to 0.01. The smoothing length was 
typically set to be slightly smaller than the mean interparti- 
cle spacing. Typical mean spacing values for 10* particles in 
units of the inner planet's semi-major axis was around 0.01, 
with our smoothing length set correspondingly to 0.01. The 
smoothing length spacing was set the same for all simula- 
tions. The timestep was set at 0.1 where the orbital period 
of the inner planet is 2tt. A timestep of 0.1 is convenient 
because it ensures accuracy while still allowing of order 10^ 
orbital periods to be completed within 24 hours for 10* parti- 
cles. Comparisons between simulations with timesteps of 0.2, 
0.1 and 0.05 yielded only insigniflcant differences in plane- 
tary migration rates, indicating that our simulations are not 
strongly sensitive to the chosen time step. All other param- 
eters are listed in table [T] The initial semi-major axes of 
the planetesimals ranges from 1.5 to 3.0. We work in units 
of the outermost planet's initial semi- major axis and with 
GM, = 1 so that the initial period of the outermost planet 
is 27r. We discuss most timescales in units of the outermost 
planet's initial rotation period. The initial inclination and 
eccentricity of the two planets were set to zero. However, 
we have found that planet eccentricities and inclinations 
rapidly settled toward zero for the planetesimal disks consid- 
ered here. All simulations were run for approximately 4000 
orbits. 



Uranus was chosen simply because in our test runs it showed 
the most interesting variation in its migration rate. Larger 
planet masses showed increasing disk disruption while lower 
masses showed greater degrees of stochastic motion. A high 
mass inner planet was initially chosen to ensure that par- 
ticles can be efflciently scattered by the outer planet and 
subsequently ejected by the inner planet. While this mecha- 
nism can function for an inner planet that has mass ratio less 
than 10~'^, the migration of the outer planet is then more 
sensitive to the distance between the planets. Our choice of 
a large inner planet reduces the dependence on the distance 
between the planets. If the distance between the two planets 
is large then it takes longer for the outer planet to scatter 
particles sufficiently that they cross the inner planet's orbit. 

The parameters of our simulations are listed in table 
[T] Three sets of simulations are run. The first ten, denoted 
Cl-ClO, include all pair-wise forces and contain a range of 
planetesimal disk masses and number of planetesimals. In 
these simulations the planetesimal mass is fixed at 10~^ of 
the stellar mass. The most massive disks simulated (simula- 
tions CI, Nl, and Vl-VlO) have a mass ratio of 10"^. This 
is a factor of a few larger than the more massive 200 Earth 
mass disks simulated by p revious studies (|Hahn fc MalhotTal 
1 19991 : IComes et al.l [iooj ). Simulations Nl-NlO are identi- 
cal to C1-C10 except forces between planetesimals are ne- 
glected. Simulations Vl-VlO compute all force pairs, have 
the same total planetesimal disk mass, but vary the number 
of planetesimals and the planetesimal mass. 

We first discuss the semi-major axis and eccentricity 
distributions during the simulations. We compare simula- 
tions with and without planetesimal self-gravity. We then 
discuss differences seen in the migration rates of the outer 
planet. Last, we consider the effects of increasing the num- 
ber of particles in a simulation while maintaining a constant 
mass disk. 



2.6 Comparison code lacking forces between 
planetesimals 

Our code can compute all interparticle forces, however most 
planet migration simulations have neglected forces between 
planetesimals. We desire a direct comparison between sim- 
ulations that take into account all inter particles forces and 
those that neglect forces between planetesimals. To do this 
we run two versions of our code that are identical except 
that interparticle forces are not computed if both particle 
masses are below that of a planet. In all other respects the 
two sets of simulations are identical. 



3 MIGRATION OF 2 PLANET SYSTEMS 

We begin our study of migration with the simplest system 
that allows migration of the outer planet into the disk; that 
of two planets just inside an initially cold (low velocity dis- 
persion) annulus of planetesimals. Our first set of simulation 
has a inner planet with mass 10~^ of that of the central star 
(similar to Jupiter) and an outer planet with a mass ratio 
of 5 X 10~^ (similar to Uranus). A low mass was initially 
chosen for the outer planet so that outer disk is not com- 
pletely disrupted during migration. A mass similar to that of 



3.1 Eccentricity and semi-major axis distributions 
during migration 

Figure [T] shows the semi-major axis vs eccentricity distribu- 
tion at different times during the CI simulation computing 
all force pairs for 10000 planetesimals. Each frame in the fig- 
ure is separated by 200 orbital times where an orbital time 
corresponds to the initial orbital period of the outer planet. 
As can be seen from the position of the outer planet in each 
frame it migrates outward during the simulation. The in- 
ner planet drifts inwards but not significantly so because it 
is more massive than the outer planet. The particles with 
semi-major axis interior to the outer planet are higher after 
the planet has passed than they were before h and, leaving a 
possi ble signature of a migrating planet (e.g., iLufkin et al.l 
120061 ). We also notice an increase in the eccentricity disper- 
sion of particles in the outer disk with time during the first 
half of the simulation. 

We can compare this simulation to an identical one, 
denoted Nl but lacking planetesimal self-gravity. The semi- 
major axis and eccentricity distributions for this simulation 
are shown in Figure (2] The increase in eccentricity disper- 
sion of particles in the outer disk seen in the self-gravitating 
disk is not present in simulation Nl. Thus the eccentric- 
ity dispersion increase in the self-gravitating outer disk is 



Figure 1. The semi-major axis (x-axis) vs eccentricity (y-axis) distribution of particles and planets different times during simulation 
CI. Each hash on the y axis is an increment of 0.2 in eccentricity starting at 0. Each hash on the x-axis is an increment of 0.5 in semi 
major axis, with the inner planet starting at 1.0 and the planetesimals in the disk taking an initial starting position between 1.5 and 3.0. 
Our migrating planet has an initial semi major axis of 1.5. This simulation has 10000 particles and computes all force pairs so that the 
disk feels self-gravity. Each frame is separated by 200 rotation periods based on the migrator's initial orbital period. The two planets are 
shown as green crosses and the planetesimals as red dots. Two scattering surfaces are seen, one associated with each planet. These are 
broad features with eccentricity and semi-major axis given a particle a planet orbit crossing periastron. In the first half of the simulation 
the eccentricities of particles in the outer disk increases due to self-gravity. Eccentricities increase for particles that are not scattered 
outward as the planet passes. 



likely caused by gravitational stirring of the planetesimals 
by themselves. 

Two dominant scattering surfaces are seen as broad 
strokes in the semi-major axis vs eccentricity distributions 
of Figure [T] and (2] one each from the inner and outer planet. 
These surfaces correspond to particles with semi-major axis 
and eccentricity that have periastron crossing a planet's or- 
bit. Similar scattering surfaces have pre viously been seen in 
simulations of mig rating planets (e.g., iGomes et al.l |2004| : 
iHahn fc Malhotral [20051 : iThommes et al.l 120081 ') and imply 
that both planets are scattering planetesimals outwards. 
However, for the outer planet to migrate outwards it must 
also scatter planetesimals inwards. The total energy and an- 
gular momentum of material scattered outwards by the in- 
ner planet must exceed that scattered outward by the outer 
planet to allow it to migrate outwards. This is one of the 
primary reasons why we required a large inner planet mass. 

When the outer planet migrates outwards it replen- 
ishes a population of moderately eccentric bodies residing 
between the planets. This allows the inner planet to con- 
tinually scatter objects. As these objects are scattered out- 



wards, they stop interacting with the outer planet allowing 
it to continue migrating outwards. Previous st udies have de- 
scrib ed migration in terms of two regimes (e.g. jGomes et ahl 
l2004h . If the outer planet migrates at increasi ng slower 
rates , its migration is described as "damped" fe.g.. llda et al.l 
i2000i ). however if it continues to migrate until it reaches 
the edge of the di sk the migration is considered "forced" 
(jGomes et al.ll2004l ). The planetesimal disks considered here 
are sufHciently massive that the planet migrates outwards 
until it is near the edge of the disk. 

A comparison between Figure [1] and [2l shows little dif- 
ference between the morphology of the scattering surfaces, 
however the location of the one associated with the outer 
planet is at a larger semi-major axis for simulation Nl than 
CI. This is because the outer planet migrates more quickly 
in the simulation lacking disk self-gravity (Nl) than that in- 
cluding self-gravity (CI). We will discuss migration rates in 
the next subsection. 

A notable difference between the simulations with and 
lacking disk self-gravity is the fraction of bodies captured 
into mean motion resonances with the outer planet. In fig- 
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Figure 2. Similar to Figure [T] except showing the simulation Nl that is identical to CI but lacking self-gravity in the disk. Each hash 
on the y axis is an increment of 0.2 in eccentricity starting at 0. Each hash on the x-axis is an increment of 0.5 in semi major axis, with 
the inner planet starting at 1.0 and the planetesimals in the disk taking an initial starting position between 1.5 and 3.0. Our migrating 
planet has an initial semi major axis of 1.5. Similar scattering surfaces are seen in this simulation. The outer planet migrates more 
quickly outwards than in simulation CI. The outer disk remains colder (has a lower eccentricity dispersion) during the first half of the 
simulation compared to the simulation with self-gravity in the disk. Populations of particles trapped in mean motion resonances with 
the planet are more prevalent in this simulation than in simulation CI. 



ure 13] we show simulations CI and Nl at 557 orbits into 
their runs (approximately 15 percent completion). We find 
that the fraction of particles captured into resonances is 
lower in the simulation with self- gravity than that without. 
While an outwardly migrating planet is capable of capturing 
bodies into mean motion resonances, the fraction of parti- 
cles captured into resonance can be reduced if the planet 
is migrating faster or i f the initial eccentricity distribu- 
tion of the disk is larger (|lda et al.ll200ol : iHahn fc Malhotral 
I2OO5I : iQuillenlliooebl ) . The lifetime of particles in resonances 
could be reduced by scattering among t he planetesimals or 
if the planet's motion is stochastic (e.g., (jHahn fc MalhotTal 
119991 )). Because the planet migration rates up to this point 
in the simulation are not significantly different in the two 
simulations, we attribute the difference in resonance popu- 
lation to either reduced capture because of a higher initial 
eccentricity distribution or because of a reduced lifetime in 
resonance. There are particles present in mean motion res- 
onances in simulation CI however the fraction of particles 
captured in resonances is small compared to that in simula- 
tion Nl lacking disk self-gravity. 



3.2 Effects of disk self-gravity and disk and 

particle mass on the planet migration rate 

In Figure 3] we show the semi-major axis of the outer 
planet as a function of time for simulations CI - CIO 
(with disk self-gravity) and for simulations Nl- NIO (lack- 
ing disk self- gravity ) . Since the simulations shown in this 
figure have planetesimals all of the same mass, the simula- 
tions with more particles have higher disk masses. We find 
that more massive disks allow the outer planet to migrate 
faster, as might be expected from scaling models predict- 
ing migration rates and previous numeric al studies (e.g., 
iHahn fc Malhotrall 19991 : [ Gomes et al.ll2004l ). In Figure [4] we 
see that there are jumps in the semi-major axes profiles with 
time, particularly for the lower mass disks with the slow- 
est migration rates. These correspond to times when the 
two planets pass through a mean motion resonance. As the 
planets are separating, the outer planet is not captured into 
resonance, rather it jumps from one side of the resonance to 
the other. For a short time the planet eccentricities increase 
and then are damped via scattering with planetesimals (a 
form of dynamical friction). The jumps in both eccentric- 
ity and semi-major axis, we see in both the self-gravitating 



7 





Figure 3. Eccentricity vs semi-major axis distribution at a time 
t = 557 rotation periods. Planets are shown with filled in circles, 
planetesimals with small dots. Scattering surfaces with perias- 
trons that are planet crossing are seen for both planets, a) For 
simulation CI including planetesimal gravitational interactions, 
b) For simulation Nl lacking planetesimal interactions. In simu- 
lation Nl particles have been captured in the 4:3 and 3:2 mean 
motion resonances at semi-major axes of 2.06 and 2.23 for the 
planet at a semi-major axis at 1.7. The outer disk in simulation 
CI is thicker and lacks as numerous a resonant population as sim- 
ulation Nl. In simulation CI either planetesimal interactions have 
scattered planetesimals out of resonances or the fraction of par- 
ticles captured into resonances is reduced because of the higher 
eccentricity distribution resulting from gravitational stirring. 



and non-self-gravitating simulations and are particularly no- 
ticeable when the disk is low mass and the migration rate 
slow. These resonant events have been previously seen in 
other numerical simulations, for example such an event is 
an important characteristic of the "Nice" model for the out- 
ward migration of the outer planets in the early solar system 
l|Tsiganis et al.ll2005l : iLevison et al.ll2008h . 

We notice differences in the migration rates in the self- 
gravitating disk compared to that lacking self-gravity. On 
average the total distance traveled by the outer planet is 
higher when self-gravity is absent, as shown in Figure [5] 
showing the total distance traveled by the outer planet af- 
ter 4000 orbits. However when we examine the total dis- 
tance traveled by the outer planet, we note that the self- 
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Figure 4. The outer planet's semi-major axis as a function of 
time, a) For simulations CI through CIO, computing all pair- 
wise forces, b) For simulations Nl through CIO, lacking disk self 
gravity. Jumps in the semi-major axis occur when the two plan- 
ets pass through resonance. The migration rate is approximately 
proportional to the disk density. Stronger migration is seen at 
the largest disk densities. Migration is not damped but neither 
is it exponentially increasing in these simulations. This suggests 
that the migration rate saturates at a level that depends on disk 
density. 



gravitating simulations have much larger scatter in this 
quantity than those lacking self gravity. The planets in the 
simulations lacking self-gravity tend to have smoother mi- 
gration excepting when the disk is very massive. 

For the most massive disks considered we note that the 
migration rate strongly increases with time for both self- 
gravitating and non-self gravitation simulations (see Fig- 
ure |4]). These may represent migration in the exponen- 
tial or "forced" reg ime predicted under some conditions by 
iGomes et al] ^2004 ). 

Following the toy model by I Gomes et al.1 (120041 ) (see 
their equations 1-3) the planet migration rate, dp is expected 
to depend on the total mass of orbit crossing planetesimals 
M{t) 



k M{t) 1 
27r Mr, 



(10) 
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where fc is a possibly migration rate dependent function that 
depends on the distribution of the planet crossing planetes- 
imal orbits, Mp is the mass of the planet, and ap is the 
planet's semi-major axis. The evolution of M{t) depends 
on the timescale, r, for the planet to scatter planetesimals 
away so that they are no longer orbit crossing (also possibly 
migration rate dependent) and the additional mass in plan- 
etesimals that becomes orbit crossing because the planet has 
moved further into the unscattered disk; 



10000 



M(t) = -M{t)/T + 2TTap\ap\a{ap), 



(11) 



where a-{a) is the surface density of the unscattered plan- 
etesimal disk. 

The sign of the parameter 



+ k^/a^a{ap)/Mp 



(12) 



determines whether the migration is damped or forced. If 
a < then M(t) decays exponentially and the planet stops 
migrating, otherwise the planet accelerates. For the most 
massive disks considered here we do see accelerations in the 
migration rate, and for the lowest density disks we see lit- 
tle migration. However for most of our simulations we see 
smooth or nearly powerlaw migration rates. For the migra- 
tion rate to fail to be exponential we require that a ~ 0. 
There could be an intermediate of disk densities, in between 
the forced and damped regimes, where fc and r can be con- 
sidered functions of the migration rate. Alternatively the 
migration could increase until a saturation level which de- 
pends on disk density. We note that the disks we consider are 
more massive than required for migration; a{r)r^ > Mp. For 
the massive disks a larger population of scattered objects is 
left after the planet passes through the disk implying that 
only a fraction of it has efficiently imparted energy and an- 
gular momentum to the planet. The planet scatters particles 
at lower efficiency during more rapid migration. This would 
lower fc and r possibly accounting for the near cancellation 
of the two terms comprising a. 

We now fix the mass of the disk and consider how the 
mass of the particles affect the migration rates. In the VI- 
VIO simulations we altered the mass of the planetesimals 
such that no matter what number of particles there where, 
the total disk mass was always a Jupiter mass. Figure [6] 
shows the distance traveled by the outer most planet and 
semi-major axis as a function of time for 3 simulations. 
While Figure |6^ shows a substantial variation in the to- 
tal distanced traveled by the planets in these simulations 
but no strong dependence on the planetesimal mass. Fig- 
ure [Bja, showing the semi-major axis profiles shows that 
the migration of the outer planet is far more stochas- 
tic when the pl anetesimal mass is high e r as explored by 
previous studies jHahn fc Malhotra|[l999l : IZhou et al.ll2002l : 
iMurrav-Clav fc Chiang||2006l ). 



4 SUMMARY AND DISCUSSION 

In this paper we have described an implementation of the 
heliocen tric democratic 2n d order symplectic integration 
scheme (|Duncan et al.|[T998l ) written in CUDA and designed 
to work on a graphics processing unit. This allows us to 
efficiently use the large number of processors available on 
the card in parallel. The code is designed to compute a 
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Figure 5. The total distance (x-axis) migrated by the outer 
planet after 4000 orbits as a function of number of planetesimals 
in the disk (y-axis). a) As measured from simulations Cl-ClO 
which compute all pair-wise forces, b) As measured from simula- 
tions Nl-NlO which neglect the self-gravity of the planetesimals. 
For these simulations the planetesimals have the same mass so 
the particle number sets the total disk mass. We find that the 
outermost planet migrates faster in proximity to a more mas- 
sive planetesimal disk at rate that is approximately proportional 
to the disk density. The planetary migration rate is on average 
~ 1.3 times higher when disk self-gravity is neglected. The lines 
show the best fit to the data points. 



large number of gravitational interactions as well as com- 
pute Keplerian advances in parallel. Because all memory 
resides on the device this parallel computation platform 
has advantages over message passage interfaces that must 
pay a penalty to share information between processors. The 
biggest drawback with our current computational platform 
is the low precision, a problem that is now being resolved 
with the recent availability of double precision capable video 
cards. 

We have used the new code to simulate planetary migra- 
tion into a planetesimal disk with 10* bodies and have com- 
puted all gravitational interactions. We have explored the 
difference between simulations of planets migrating into self- 
gravitating planetesimal disks and those lacking computed 
interactions between planetesimals. We find that the fraction 
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massive disks. However for the most of the disk densities 
considered no exponential increase in migration rate is seen 
implying that the migration rate saturates at a level approx- 
imately proportional to the disk density. 
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Figure 6. a) The total distance traveled by the outer planet 
in 4000 orbits as a function of the number of planetesimals in 
the disk for simulations Vl-VlO. For these simulations the total 
disk mass is held fixed so that the the planetesimal number is 
inversely proportional to their mass. The total distance traveled 
by the outer planet is not strongly dependent on the mass of the 
planetesimals. b) The semi-major axis as a function of time in 
orbits for simulations VI, V6, and VfO, with the same disk mass 
but different masses and numbers of planetesimals. The three sets 
of data points have semi-major axes offset by 0.3 so the three tra- 
jectories can be seen together on the same plot. The key lists the 
simulations and planetesimal mass ratios for those simulations. 
The planet motion is more stochastic when the planetesimals are 
more massive. 



of particles present in mean motion resonances is reduced 
in the simulations with self-gravitating disks compared to 
those lacking planetesimal interactions. We attribute this to 
a reduction in the resonance capture rate due to the larger 
eccentricity dispersion caused by gravitational stirring in the 
self-gravitating planetesimal disks. This suggests that disks 
in proximity to migr ating planets coul d be featureless (as is 
true for Fomalhaut. iKalas et al.ll2005h . A dust disk lacking 
resonant structure could still host migrating planets. 

We find that planet migration is smoother but some- 
what faster in the simulations lacking self-gravity. Migra- 
tion rates in both cases are approximately linearly depen- 
dent on planetesimal disk density. The planet can undergo 
rapid changes in migration in the simulations with the most 
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Table 1. Simulations 
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Here Mi, M2 are the masses of the two planets, and m is the planetesimal mass in units 
of the central star mass. The parameters ai and 02 are the initial planet semi-major axes. 
The parameter N is the total number of planetesimals, and is the total planetesimal 
disk mass in units of the stellar mass. The planetesimals initially extend between 1.5 and 
3.0 in semi-major axis. Simulations Cl-ClO included all pair-wise forces all have the same 
planetesimal masses but a range of total disk masses. Simulations Nl-NlO arc identical to the 
C series but neglect interactions between the planetesimals. Simulations Vl-VlO compute all 
pair-wise forces and vary planetesimal mass in order to maintain constant mass disk with a 
fixed number of particles. Simulations were run for approximately 4000 orbits. 



